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1 Introduction 

We consider the problem of finding a matrix U 6 3v nxn such that 



U T (T-SX)U = A- IX, 



is diagonal, or equivalently 



U T SU = I and U t TU = A, 



( 1 ) 



where 
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ft- 1 




7n - 1 


ft — 1 ^ n 
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and 5 is assumed to be positive definite. This generalized eigenvalue problem has 
two special cases that are of interest in themselves. They are: 

1. 5 = 7, the ordinary tridiagonal eigenproblem. 

2. 5=7 and a* = 0, the bidiagonal singular value problem (bsvp), by perfect 
shuffle of the Jordan matrix 



0 B t ‘ 
B 0 



with B upper bidiagonal [16]. 

There are two phases to the divide and conquer algorithm, the divide (or split) 
phase, and the conquer (or consolidate) phase. We shall address these in order. 



2 The algorithm 

2.1 The divide phase 

Denote by e,- the ?th axis vector where the dimension will be clear from the context. 
Let s, 1 < s < n, be an integer, the split index , and consider the following block 
forms: 





Ti 


ftj-ift-i 




T = 


A-ieTLi 


ft* 


ft c T 






eift 


T 2 




Si 


e 5-i 7.^-1 




5 = 


7 s — 1 e J'— 1 


ft 

e l7a 


is 

s 2 



Note that s = n is possible; then T 2 , 52, and ej are empty [9, 10]. Suppose we 
solve the subproblems 



UZ(T k -S k \)U k = A k -I\ {k = 1,2). 



( 2 ) 
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The form of the subproblems is preserved. In particular, the matrices Sk are 
positive definite and, if T lias a zero diagonal, so do the matrices Tk . Let 



Then 



0 



U i 

1 

U 2 



U T {T -SX)U = 

' UT(T 1 -S 1 X)Ui £/ 1 T e J _ 1 (A_ 1 - 7 ,_ 1 A) 

(A-i ~ 7 s-iX)eJ_ 1 U 1 a s -6 s A (ft - 7,A)e|’t/ 2 

t/ 2 T ei (A - 7» A) t/ 2 T (T 2 - 5 2 A)f/ 2 



2.2 The conquer phase 

The conquer phase consists of solving the subproblems (2) from the divide phase, 
consolidating the solutions, and finally, solving the consolidated problem. Let 

Ui = u I Q,-i, u 2 = f/Jej, 



where the Uk are solutions to (2). Then 

Ai - IX 

U T (T- SX)U = ' 



ui(ft-i - 7 s-iA) 

a s - 6 S X (A— 7«A)uT 
Ao - /A 



(A-i - 7*-i^)«r 

« 2 (ft - 7s A) 

The right side is the sum of a diagonal and a Swiss cross : 



U T ( T - SX) U = 



This can be permuted to an arrow matrix by a permutation similarity transfor- 
mation with P s = [ei,e 2 , ...,e,_i,e 5+ i, ..., en,e 5 ]. Thus 



X 
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X 




X 
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X X X X X 
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j4(A) := PjU T (T- SX)UP 3 
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with 
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n = 



B = 






a/ 



ui 

ll 2 J ’ 

C = 



1,-1 1 



1,1 



Since S and 



7 Cu 
n T C 7 

from the former. Its Cholesky decomposition is 



are congruent the latter inherits positive definiteness 



= R t R, 



'I Cxi' 




7 




‘ 7 


Cu 


u T C 7 




u T C p 






P 



with — 7 — ir C“u > 0 the Scliur complement in 5 of 






Now 



R~ l = 



I —Cu/p 

l /p 

and a second congruence transformation with gives 



= R~ 



- T A{X)R- 


i— T 


D 




u t B 


D 


w 


T 


U) 



with 



7T 1 - IX 
-IX =: A - XI 



(B - DC) u 

P 

o, — u T (2 B — DC) Cu 
P 2 



We have reduced the conquer step to the problem of solving an ordinary eigen- 
problem for a symmetric arrow matrix. If V is an orthogonal matrix with 

AV = VA 

and A diagonal, then (1) holds with 

U = U P. R~ 1 V 



U\ -Ui un,-i/p 

1/P 

U-> -UnW-ilJp 



V. 
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It is useful that v* = UkUk can be computed in 0(7i) time by solving S\V\ — 
e 5 -i and 52 v 2 = ei using the Cholesky factorization S\ = L\Lj and the reverse 
Cholesky factorization 5 2 = In the case that only the eigenvalues are 

wanted it is only necessary to compute the first and last rows of the [/-matrices 
which constitutes a further savings. 

In summary, the conquer phase proceeds by consolidating the subproblems and 
building a full eigenproblem for an arrow matrix. 



3 Solving the eigenproblem for the arrow 

In this section we consider the solution of the eigenproblem for a real symmetric 
arrow matrix 



A = 



D b ‘ 
b T 7 _ 



where A 6 3? nxn is symmetric, D = diag(a), a = [aq, ...,a n _i] T , ct\ > a 2 > ... > 
a n _i, and b = [f}\ , ...,/? n _i] T > 0. When A arises from the BSVD then a is odd 
and b is even, that is a + Ja = 0 and b = J b, with J the counter-identity , the 
identity matrix with its columns reversed, and 7 = 0. 

If any (3j — 0 then it is possible to set A j — a j and deflate the matrix since 
ej is clearly an eigenvector [28]. We shall call this /3-deflation and note that if 
(3j < /o/^||b|| where tolp is a small tolerance then a numerical deflation occurs. 
We derive a precise value for tolp in section 4.4. 

A second type of deflation occurs if applying a 2 x 2 rotation similarity trans- 
formation in the ( j, j + l)-plane that takes fij to zero introduces a sufficiently small 
element in the (j, j-f 1) position of the matrix. This will be called a co?n6(>-deflation 
(see [15]). At each consolidation step we perform a sweep to check for //-deflations 
followed by a sweep to check for ccnuto-deflations. The a>?n6o-deflation can be 
arranged so that the ordering of the c*j is preserved whenever one occurs. After 

deflation the new /?J +1 := yjfij -f /?? +1 > /?y +1 and hence no further /3-deflation 

can occur. The amitodeflations can be disposed of with a single pass by backing 
up a single element whenever one occurs. Note that deflation is backward stable 
since it results in small backward errors in A. Deflation for the BSVD is more 
delicate involving a simultaneous sweep from both ends of the matrix. Care must 
be exercised at the center of the matrix. 

After deflation the resulting matrix can be taken to have all fij > 0 and the 
elements of the arrow shaft distinct and ordered, that is a\ > c*2 > ... > a n -i* 
An arrow matrix of this form will be called ordered and reduced. Henceforth, we 
shall assume A is of this form. 

The block Gauss factorization of A — XI is 
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A -XI = 

where /, the spectral function of A, is given by 



I 0 ' 


' D- XI 


b 


b T (D - A/)" 1 1 _ 


0 T 


-m . 



n - 1 



/(A) = A - 7 + £ -A 

' rv • — 



j =i a ;' 

This is a rational Pick function with a pole at infinity [1]. The most general form 
of a rational Pick function is 



n- 1 

/(A) = 6X- 7 + £ 
j = l 




(5 > 0 . 



(3) 



In relation to the various divide and conquer schemes, the case <5 > 0 corresponds 
with extension , 6 = 0 with modification , and 6 = 7 = 0 with restriction [7]. 

Inspection of the graph of the spectral function reveals that the elements of 
the shaft interlace the eigenvalues 



Ai > cm > A 2 > ... > a n _! > X n . (4) 

Moreover, in the present case, the derivative of the spectral function is bounded 
below by one so that its zeros are, in a certain sense, well determined. 



3.1 The zero finder 

The fundamental problem in finding the eigenvalues of an arrow is that of providing 
a stable and efficient method for finding the zeros of the spectral function. We 
now examine this problem in some detail. 

The zero finding algorithm we present is globally convergent in the sense that 
the iteration will converge to the unique zero of / in (cu ,a*_i) from any starting 
value in the closed interval [o/,., a*_i], where we put a 0 = +00 and a n = — 00 . 
The zero finder converges inonotonically at a cubic rate and applies to a general 
Pick function as given in formula (3). 



3.2 Interior intervals 



The iterative procedure for finding the unique zero of / in one of the interior 
intervals ( 0 *, n*-i) proceeds as follows. Let xo, a* < xo < a*-i, be an initial 
approximation to A* . If Xj is known choose 



<j) 2 (x) - a + 



U>Q 



+ 



^1 



so that 



a Jfc-i - * 



a k - x 
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^ i \r j ) = f^(x j ), i= 0,1,2. 
Thus o’, u>o, and uq must satisfy 



(5) 



1 (or/t -1 - Xj) 1 (at - Xj) 1 




G 




/(Xj) 


0 ((\k-\-Xj) 2 ( a k - Xj) - 




Uq 


= 


/'(Xj) 


_ 0 (a k -i - Xj)~ 3 (a* - Xj)~ 3 _ 




. ^1 




L /"(*>) J 



and we find 



a 



Uq 



Ui 



_ / v v — ^ Pi O* — Or — 1 Of — Of Jt 

3r j -( 7 + a,_ 1 +a,)+ V — — L 

(V; — Xi a, — Xj c\i — Xj 

ija-i.i- ' ; 1 1 1 1 



Pi~i + 

Pi + 



( a t - 1 - Xj) 3 



— O’t 

(X; - a,) 3 



i+ , * v" 1 v* i • 



p’i «. - a/;. 



Ofc-l - Ofc 



1+ H t _ r 

i*h?I,k( ai ai Xj 



Since cjo > 0 and uq > 0 it follows that <j) 3 is a Pick function. Thus <j>j has a 
unique zero £j +1 G (a*, Also 



> #L, > o , uq > $ > 0. 

The error function 



f(x) - <j>(x) - x - (7 + a) + ^ 

t ?k-l,k 




+ 



ffi-1 ~ 

«A— 1 - X 



t 4-^1 

a k - x 



has n zeros, counting multiplicities. There are ?/ — 3 zeros exterior to (cv*, or* — i ) and 
three more at x.j. Thus t he error function crosses zero exactly once in the interval 
Hence Xj + \ lies between xj and A*, and the iteration is monotonically 
convergent from any starting guess a: 0 G [ct * , o * _ i ] as claimed. The cubic rate of 
convergence follows from (5). 

Successive iterates can be found by solving quadratic equations. Rather than 
solve <t>j(x) = 0 for £j+i it is better to solve 



M*J - A) = 0 

for the increment A = Xj — xj + Some rearrangement using (5) reduces this to 

a A 2 + /iA — / = 0, (6) 



with 
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a = 



p 



(ojfc-1 - 

= f'( x j ) — 



- n k)’ 

1 

+ ' 



( 7 ) 



1 



/(*>)• 



C*k ~ Xj J 

When shifts of the origin to the nearest pole [15] are used then one of a*-i or at 
is zero. The computation of [3 = P(x.j) should account for the fact that it has only 
simple poles at ak-i and a*. 

If we start at the midpoint of the interval, xq = -fa*)/ 2, then we always 

have fi — f3(xj) > f(xj) > 1. This can be seen by noting that /3(x 0 ) = f'(x 0 ) 
and that when xq > A*. then for all of the succeeding iterates f(xj) > 0, by 



monotonicity, and 



+ 



is negative. If xq < A* a similar argument 



~ x j 

applies. It follows that the increment can always be computed stably as 



A = 



VIP 



i + 



/l 4- 

V + Tt 



(9) 



3.3 Exterior intervals 

The treatment of the two exterior intervals is geometrically the same as above. 
Again, the approximating function has poles at the endpoints and the residues at 
these poles, and the constant term, are chosen to satisfy (5). We present the case 
for the interval (ai,oo), the case for the other exterior interval being similar. Now 



with 



0j( x) = u Q x - a + 






Cl \ — X 



n- 1 



U/'Q 



= ■ + £ , l>} P"~"‘ 2 1 

fit ( J *j ” «i)“ Xj - C*i 



71 - 1 



- = * + 2 *. 

The inequalities are strict unless n — 2. Again we find (6) where now 

Mi— 1 0? 






l 4- e, tn-ft, 

1 ' Z-^t = 2 (X 3 -O t ) : X^-a, 

Xj - O] 



P = 

nr* 



/(*j) 



These are limiting cases of (7) and (8) ( introduce another pole c*o > c*i and let 
a 0 — * +oo). If x*o > Aj then /(A) > 0 so ft > f > 1 and A is again computed 
stably using (9). We obtain global monotone cubic convergence as before. 
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Contrary to the algorithms of [11, 12, 10] our algorithm is well-defined when 
starting at the endpoints of the intervals. The algorithm of [23] can start at the 
endpoints but has only quadratic convergence. 

To guarantee that xq > X\ we take xq to be the iterate in (ai,-foo) from Too. 
As x 0 —> +oo the approximate Pick function tends to 



G i — X 

Our actual starting guess is the zero of (10) in (gi,+oo): 



( 10 ) 






*0 = < 



O] + 



When shifts are used we have c\\ = 0. 



+ IIHl" i i > oi, 



i 7<^i- 



+m j 



3.4 Orthogonality of the eigenvectors 

It is essential that the computed eigenvectors of the arrow matrix be numerically 
orthogonal. As a point of entry into the further analysis of the algorithm we now 
examine the orthogonality of the eigenvectors following [15]. 

Consider the divided difference 



/(A./0 



/O) — /(/') _ | y- @j_ 

X-fi ~ - A )( a ; - /') 

l + b r (D-A/)- 1 (D-///)- 1 b. 



Note that /< = A gives /'(A) = ] + ||(Z7 — A/) 'b||r'. If /(A) = 0 then 



(11) 



r(A) = 



07 

A — ti j 




1 

I 


l 




l 



is an eigenvector of the arrow matrix A 
value A, and 



D l» 



7 . 



associated with the eigen- 



«(A) = 



-■(A) 

\/m 



is the normalized eigenvector whose last element is positive. The ordering of A 
implies that its matrix of eigenvectors can be taken positive below and on the 
diagonal, and negative above. 
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Let /( A 0 ) = f(iio) = 0 with A 0 ^ //o- Thus A 0 and / iq are distinct eigenvalues 
of A. The eigenvectors u(Ao) and v(no) are orthononnal: 

h(A 0 ) 7 u(/i 0 ) = /(A 0 , /<o) = 0. 

Let A and // be approximate eigenvalues in the sense that 



A- Aq 

C\j — Aq 




y _ /< ~ /'U 

3 - /'(/ ’ 




( 12 ) 



Here 6 > 0 is hopefully, but not necessarily, close to the machine unit c. Note that 
(12) is equivalent, with 



- A 



— 1 T &j , 



~ /' _ 



— = 1 + tf. 

— I In J 



L\j — A 0 ' f\j - no 

These conditions imply that the approximate eigenvectors u(A) and u(ji) are 
nearly orthogonal. For we have 



vWtf'O'MA ) 7 «(/i) 



Since 



/('V I 1 ) ~ /(AuW^o) 

n - 1 

E 






J = 1 
n- 1 

£ 

j=i 



(<>j - A)(Oj - n) 

$ 

(<‘j - A)(ttj - //) 



- A)(ftj - /<) \ 
( n j - Ao)(«j -tto)J 



(Sj + S'j + SjS'j). 



Ify + 6j + 6j6j I £ 



26 



+ 



8 - 



1 + 6 (l + d) 2 



< 2 6 



then 

/(W^'UA) 7 ' «(//) = 261 7 (D - A/)-'0(D - y/7)- ] b 
with |© | < /. Thus 



v//'(A)/'(/0MA) T t,(,,)| = 26||(D- A/)- 1 l,|| 2 ||(D-/i/)- , b|| 2 , 

and so 



|h(A) t i/(//)| < 26. 

Condition (12) is stringent. If we let /3/ ; — 0 then it is easy to show that A 
can have an eigenvalue A 0 = Ao(,4) = (u- + 0(/ijj:); (12) then reqiiires that the 
approximate eigenvalue A satisfies a bound 
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|A- A 0 |<O(*ft 2 ), 

which is difficult if ft./||l>|| is only somewhat larger than machine precision, say 
e 3 /4 r p vvo techniques are used to attempt to satisfy (12) - shifts of the origin [15], 
and simulated extended precision (sep) arithmetic [26, 14]. Condition (12) means 
that 



|A - A 0 | < 6 min{A 0 - <U--i - A 0 ). 

When shifts are used it means that A is nearly //( Ao). 



4 Numerical stability of the algorithm 

We now give a partial analysis of the .stability of this approach to the eigenproblem 
for the symmetric arrow matrix. Observe that 

, m = £ w n;.i(*--v) 
iW n ;r,'(A-..j)' 

The following inverse eigenvalue probh w [6] is important.: given {cv ; ) and {Aj} 
satisfying (4), find {ft} and 7 so that A(/l) = { A j } . This problem is simply solved 
by computing the residues of the partial fraction decomposition of /. In particular 



ft J 



1 



lim (n*. - A) 

A — 



HA) 

'/(A) 



_ ~ ah 

n n — I 

j = 1 ;=i 



For fixed { cvy } , the elements of the arrow head, {ft} and 7, are explicitly known 
functions of the eigenvalues. 

Now let {Aj} be a set of appvorimat c eigenvalues of A satisfying (4). Then 



Hi 



rij‘=i( f H- - Aj) 

Ylj?k( c 'k - <v) 



(Hk > 0), 



(13) 



n n — 1 

7 = (n) 

; = 1 3 - 1 

define a modified matrix A with A(.4) = {Aj). To obtain a backward error analysis 
for the complete eigenvalue problem we bound the differences ft — ft. and 7 — 7. 
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4.1 Error analysis for the Dongarra-Sorensen condition 

We give an error analysis using (In' Doiigarra-Sorensen condition 

= M<«, (15) 

t\ k Aj 

where 6 = 0(c) is of the order of t lie machine unit, simplifying that in [6]. 
Rearrangement of (15) gives 

Aj - a* = (Aj - n*)(l + 

It follows that 

n { n 

~ fy.*) ~ 1 + XI 

j= 1 \ j=l 

and 




A - ‘h- ^ 1 + 6 j,k 
with the Sj k and 6" k at most, only slightly larger than the Thus 



where 6" — O(c) is only slightly larger than 
Now (M) become* 




7 — 7 + X^^ J n A(j ) )^j,k(j ) 

J = i 

with a*(j) one of the poles of /. d Inis 

71 

b ” 7 1 £ ^ XI l A -? “ a ^(j)l- 

j = i 

To minimize this hound we choose n^j } to he a pole of / closest to Aj. Clearly, 
a *(i) = and n A(n) = <»„_!, so 

Id _ d| < ^ A 1 - o ! ) + I Aj - C»Jt(j)| + (<>n-l - A n )^ • 

For 1 < j < n a closest pole to Aj is either Oj or Oj_i. The distance 



|Aj ~ n k-[j )l = niin {Aj - aj^vj.j - A j } 
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is maximized when Xj is the midpoint, of the interval (cij,cij_ i), and the value of 
the maximum is (aj + Oj_i)/2. Thus 



n- 1 



It — tI < 6 f (Ai — O 1 ) + — y (nj_i - (Ij) + - A n ) 



j=-' 



< 6{X } -Xn) < '26\\A\\->. 



In summary, the Dongarra-Sorcnsen condition implies small relative errors in 
each Pk and a small absolute error in y . For the BSVD this implies small element- 
wise relative errors since the condition — - — 0 is enforced by Xj + A ri+ i_j = 0 
(only half of the eigenvalues are actually computed, the rest follow from this con- 
dition). 



4.2 Rounding error analysis of the computation of /(A) 

The choice of a termination criterion depends on a careful rounding error anal- 
ysis of the particular manner in which we compute /(A). Let {a,}, {Pi}, and 7 
be floating point numbers. We rcpn^mt X as the ordered pair of floating point 
numbers (<r, // ) where the shift a is a pole closest to A, and A := cr -f //. For the 
exterior intervals we have a — n 1 or a — n„_j. For the interior intervals cr can be 
determined bv evaluating f at tin* midpoint and checking the sign. We compute 
/(A) as 



1 

foil 1 ) = T —r 1 1 - (/' ~ */')■ 

U ~ '' 



with the standard operation precedence rules, when 1 



cx'j — c\j — a and ~ — a. 

We use Wilkinson’s notation: fl(u: + 1 /) = (./* * //)( 1 + 6 ) with \6\ < t/(\ + c) 
and c — the machine unit. More generally, < denotes numbers not essentially 
larger than 2“* [27] and the rounding errors 6 satisfy |<N| < c . 

We define 



//(a . „ A ) := //(o' -;/) = //((«, - a)-//). 



If cr = 0 *. then 



fl(c\ k - A) = -// = a/,. - A, 
with no rounding error. For j ^ /•, 
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/*(«> - A) = (i\j - A) I + 



6 + , 



and since a* is a pole closest to A then 
computed with small relative errors; 



o , — A 



~ in 
o j — A 

< 2 . Thus all terms c* ; — A are 



J l{dj — A) — (fij — A)( 1 + ), |<5j | < e. 



(16) 



When computing /(A) = fA/i) "e add the term A — *) = (A — a) — (7 — cr) last. 
A routine error analysis using (16) and 

lA-^ISI/WI + f 

to eliminate the term |A — - ; | from the error hound gives 

( " _1 0 ? \ 

l/W)) - /(A)| < c I :i|/(A)| + - A| + (» + 5) £ 

which implies 



l/f(/(^))l < 0 + )|./ ('Ml ■+■ 1 ( k “ A| + (?< + 5) ^ 



p] 



la; — A I 

J =1 1 1 \ 



(17) 



4.3 Termination 

Our go;tl is lo choose a terminal ion criterion so that we stop when A is as close 
to the true eigenvalue A/, as po^.»ih|. Let // = At in ( 11 ) with /(At) = 0. Now 
o* < A h < o;-_j. Also h’t o|. • A < o/,_|. Then the terms a j — A and a j — At 
have the same sign and 



I A - At | < 



I ./ ( A ) | 



+ r;: 



(18) 



j - i a(|uj-Aa | 



To obtain an upper hound for |A — A/-| we need an upper bound for |/(A)| and a 
lower bound for the denominator. For the latter we have 



n — I 



J = 1 



- 1 d; 

~ t joj-Aj 

\ Cy j — A 1 1 r i ^ — A | ~~ ‘ inaxj |«j - Afc| 






> i + 



(19) 



Let us determine how small |/(A)| is when A is the rounded representation of 
Ajt. This is 



Ceneinlized Divide and Conquer 



and we have 



Thus 



A := a + fl(ftk) = o’ + /a- ( 1 4- <*) 
“ A* -f = Ay- + (A* — rr)^ 



|a — Ay. | — nun |n — Ay ; |. 
j 



n — 1 



|/(A)| = |A-A,| 1 + 53 



a i 



J = l ( n j ~ A)('»j - A *) 



n- 1 



— |(tf — A;,. )< L '| [ 1 + 



ft? 



JTi ( n J - A )(“j “ At) 



< ( I |(T — A;,. | + ^ 






n I'v - a i 



From (17), 



N - I j'J 

l//(/(A))l < < ( k - At I + |A - n\ + („ + (i) 53 — ^ 



l a J - *1 



Since A*. — cr = (A - cr)/(\ -f b) then 



l/^(i(A))| < t | 2|A - rr| -f (// -f 6) ^ - 



J = i 



;-^l 



We terminate and set Ay. : = A when 



/ " - 1 yj 

l//(/(A))| < * f 2 |A - <r| + (,» + (i) ^ — f- A , 
Inequality (17) also holds if /(A) anil //(/(A)) are interchanged. Thus 

( n_l A- 

l/(Afc)l < < A | A;.. - it| + (;(</ + 17) 53 ~r . 

\ >=iK“At| 



From (18) and (19) 



- Ai . i +(:<!, + 17)^;:;^ 



I A t - At | < i max |A,, - i\j \ 
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Since |cr — A* | < \cr — A/.. | 4- | X t , — Ay. | and |<r — A/. | < maxj |a ; — A*| the computed 
eigenvalues satisfy 



| A/,. - A/.. | < (dn + I7)t max |c»j - Ajt | 

j 

< (i(u + G)< 



( 21 ) 



4.4 Error analysis for the Gu-Eisenstat condition 

From j — 7 = i ( Aj — Aj ) and (21) we find 

h - 7 I < + ())<||d|| 2 . 

We have noted that the Ihnigarra-Sorensen condition (15) is stringent. It is 
natural to ask for small ab.suhil, « nors in tin* ,i f: . if we replace 6j j. by Sj } k/fi k i n 
the analysis in section 1.1 we iind that 



— >h: I 4- 



i "~ 1 a" \ i ri_1 

-Xf 

3 - 1 J j=1 



and 



M, -d,|< n < *(i 4-0(0), 



are implied by the (iu-Kist n.sltil unnhhun 



A, - A, , 

-d/,. — = < d. 

'»/ - A 



We must bound b. 
From (20) 



u - 1 






I A * — Ajt | I 1 + 5 ) — 1Ui ( I A/. ~~ a \ + 5Z "i ” 

V -A/,)(n 7 - A,jy ^ l a > - 

with m = 3(?/ 4- 6) Using 

I A — « T j < | A — Ay. | 4- | A/. — a | 
and the Gu-Eisenstat inequality. 

i/-' 



A* | 



1 



K - a| 



< 



(a j - A)(oj - Ay. ) 



4- 



|A, - A| 



{(Xj - A )( ci j - Ajt) * 



we get 



(iCMU'raliml Divide and Conquer 
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|A* - A* | < 



|Ac- - <t| + 



r;: 



- 1 






r;: 



1 [Oq-Aj )( n j - Aj. )] ] t'* 

~\ 



1 luj-A* K^j-Afc) / 



where t has been increased to c/( i - im ). 
By Cauchy’s inequality, 






|A*-A*| < me 



| A/; - (T| + 



\ 






f 'll ) 1/_ i 

(rij — Aj. )(« _, — Ajl ) / f 

< UU ( |A/; - CT\ -f |(nj - A/; )(Uj - A/; ) 



for every j. The arit liniei ic-geomet ric mean i 1 1 1 q u : 1 1 i i \ and the triangle inequality 
yield 



|Afc — A/ : | < im ^|A/ ; — cr\ -f ^ | nj — A/,. | -f — | A /■ — A; : | 

Thus 



Al . , , A ,IN\ 

2 — ~ J h* ” - /m U^ /; ~ ^1 + hi “ ^ /; l”^~ ) 

. I A/.- - rr| /?j 



, 1 1 hi “ 

— /j /( j| 

h 

< ~im< [|l i|| ---■■^-—■1 



|A/,-^|||b|| 



If me||b|| < )3j for all j, then 



hi,/* I < ■'!///( ||1>||. 



and consequently 



\>h- “ h.| < G//(// + ())< ||b||. 

Thus tolp is me. If /i*. < 3(?/ 4- G)(|]/d| wa replace ,i/. by zero and accept £is 
an eigenvalue with normalized eigenvector <»/... 

The computed eigenvectors of .4 are tala n to lx* those of the nearby matrix A. 
Because of (13) and (16) they are computed to high relative precision elementwise 
and hence are numerically orthogonal [20]. 
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